INVERSE LAX-WENDROFF METHOD FOR BOUNDARY CONDITIONS 

OF BOLTZMANN TYPE MODELS 



FRANCIS FILBET AND CHANG YANG 



Abstract. In this paper we present a new algorithm based on a Cartesian mesh for the 
numerical approximation of kinetic models on complex geometry boundary. Due to the high 
dimensional property, numerical algorithms based on unstructured meshes for a complex 
geometry are not appropriate. Here we propose to develop an inverse Lax-Wendroff pro- 
Q f cedure, which was recently introduced for conservation laws [21], to the kinetic equations. 

D ■ Applications in ID x 3D and 2D x 3D of this algorithm for Boltzmann type operators (BGK, 

^0 ' ES-BGK models) are then presented and numerical results illustrate the accuracy properties 

of this algorithm. 
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1. Introduction 

We are interested in the numerical approximation of solutions to kinetic equations set 
in a complex geometry with different type of boundary conditions. Unfortunately, classical 
structured or unstructured meshes are not appropriate due to the high dimensional property 
of kinetic problem. In contrast, the Cartesian mesh makes the numerical method efficient 
and easy to implement. The difficulty is that obviously grid points are usually not located 
on the physical boundary when using a Cartesian mesh, thus a suitable numerical method 
to capture the boundary condition on the complex geometry is required. Several numerical 
methods based on Cartesian mesh have been developed in computational fluid dynamics in 
last decade. Among these methods, the immersed boundary method (IBM), first introduced 
by Peskin [17] for the study of biological fluid mechanics problems, has attracted considerable 
attention because of its use of regular Cartesian grid and great simplification of tedious 
grid generation task. The basic idea of immersed boundary method is that the effect of 
the immersed boundary on the surrounding fluid is represented through the introduction of 
forcing terms in the momentum equations. In conservation laws, two major classes immersed 
boundary like methods can be distinguished on different discretization types. The first class is 
Cartesian cut-cell method [12], which is based on a finite volume method. This conceptually 
simple approach "cuts" solid bodies out of a background Cartesian mesh. Thus we have 
several polygons (cut-cells) along the boundary. Then the numerical flux at the boundary 
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of these cut-cells are imposed by using the real boundary conditions. This method satisfies 
well the conservation laws, however to determine the polygons is still a delicate issue. The 
second class is based on finite difference method. To achieve a high order interior scheme, 
several ghost points behind the boundary are added. For instance for solving hyperbolic 
conservations laws, an inverse Lax-Wendroff type procedure is used to impose some artificial 
values on the ghost points [21]. The idea of the inverse Lax-Wendroff procedure (ILW) 
is to use successively the partial differential equation to write the normal derivatives at the 
inflow boundary in terms of the tangential and time derivatives of given boundary conditions. 
From these normal derivatives, we can obtain accurate values of ghost points using a Taylor 
expansion of the solution from a point located on the boundary. 

The goal of this paper is to extend the inverse Lax-Wendroff procedure to kinetic equations 
together with an efficient time discretization technique [5, 6] for problems where boundary 
conditions play a significant role in the long time asymptotic behavior of the solution. In 
particular, for low speed and low Knudsen flows for which DSMC methods are unsuitable. 
Therefore, the main issue relies on that the inflow is no longer a given function, while it 
is determined by the outflow. For this, we proceed in three steps: we first compute the 
outflow at the ghost points. To maintain high order accuracy and to prevent oscillations 
caused by shocks, we use a weighted essentially non-oscillatory (WENO) type extrapolation 
to approximate the ghost points by using the values of interior mesh points. In the same 
time, we can extrapolate the outflow located at the boundary associated with ghost points. 
Then, we compute the inflow at the boundary by using the outflow obtained in the first step 
and Maxwell's boundary conditions. Finally, we perform the inverse Lax-Wendroff procedure 
to approximate the inflow on the ghost points, where the key point is to replace the normal 
derivatives by a reformulation of the original kinetic equation. 

For simplicity, we only consider simple collision operator as we adapt the ellipsoidal statis- 
tics BGK or ES-BGK model introduced by Holway [9]. This model gives the correct transport 
coefficients for Navier-Stokes approximation, so that Boltzmann or ES-BGK simulations are 
expected to give the same results for dense gases. Let us emphasize that Direct Simulation 
Monte-Carlo methods (DSMC) have been performed to the ES-BGK model in complex ge- 
ometry. However DSMC approach is not computationally efficient for nonstationary or low 
Mach number flows due to the requirement to perform large amounts of data simpling in 
order to reduce the statistical noise. In contrast, F. Filbet & S. Jin recently proposed a 
deterministic asymptotic preserving scheme for the ES-BGK model, where the entire equa- 
tion can be solved explicitly and it can capture the macroscopic fluid dynamic limit even if 
the small scale determined by the Knudsen number is not numerically resolved [7]. We will 
use this scheme to solve ES-BGK model while on the boundary the inverse Lax-Wendroff 
procedure will be applied. 

The outline of the paper is as follows. In Section 2 we describe precisely the inverse Lax- 
Wendroff procedure to Maxwell's boundary condition in ID and 2D space dimension. Then 
in Section 3 we present the ES-BGK model and the application of inverse Lax-Wendroff 
procedure to this model. In Section 4 a various numerical examples are provided in ID x 3D 
and 2D x 3D to demonstrate the interest and the efficiency of our method in term of accuracy 
and complexity. Finally a conclusion and some perspectives are given in Section 5. 

2. Numerical method to Maxwell's boundary conditions 
The fundamental kinetic equation for rarefied gas is the Boltzmann equation 

(2-1) ^ + vV x / = iQ(/), 

at e 

which governs the evolution of the density /(t,x,v) of monoatomic particles in the phase, 
where x e Q x c M dx , v e M 3 . The collision operator is either given by the full Boltzmann 
operator 
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(2.2) G(/)(v) = jf s ^fl(|v-v*|,cos0) (/!/' - /*/) dadv* 

or by a simplified model as the BGK or ES-BKG operator (see the next section). Boltzmann's 
type collision operator has the fundamental properties of conserving mass, momentum and 
energy: at the formal level 

f Q(f) 0(v) dw = 0, 0(v) = l,v, |v| 2 . 
Jr 3 

Moreover, the equilibrium is the local Maxwellian distribution namely: 



/ |u(t,x)-v[ : 
\ 2T(t,x) 



where p, u, T are the density, macroscopic velocity and the temperature of the gas, defined 
by 

p0, x ) = / , /(*>x,v)dv, 
Jr 3 

(2.3) ^ u(t,x) = -^- y ^ 3 v/(t,x,v)dv, 

T(t,x) = — ^- £ |u(t,x)-v| 2 /(t,x,v) a !v. 
3p(i,x) .VR 3 

In order to define completely the mathematical problem (2.1), suitable boundary conditions 
on <9f2 x should be appled. Here we consider wall type boundary conditions introduced by 
Maxwell [15], which is assumed that the fraction (1 - a) of the emerging particles has been 
reflected elastically at the wall, whereas the remaining fraction a is thermalized and leaves the 
wall in a Maxwellian distribution. The parameter a is called accommodation coefficient [4]. 

More precisely, for x e 9f2 x the smooth boundary <9f2 x is assumed to have a unit inward 
normal n(x) and for v • n(x) > 0, we assume that at the solid boundary a fraction a of 
particles is absorbed by the wall and then re-emitted with the velocities corresponding to 
those in a still gas at the temperature of the solid wall, while the remaining portion (1 - a) 
is perfectly reflected. This is equivalent to impose for the ingoing velocities 

(2.4) /(t,x,v) = (l-a)ft[/(t,x,v)] + aM[f(t,x,v)], xe^, v n(x) > 0, 
with < a < 1 and 

ft[/(i,x,v)] = /(t,x,v-2(vn(x))n(x)), 
M[f(t,x,v)] = /z(t,x)/™(v). 
By denoting T w the temperature of the solid boundary, f w is given by 



(2.5) 



(2.6) / w (v) := exp - 



^ ± u, 

and the value of /i(i,x) is determined by mass conservation at the surface of the wall for any 
t e M + and x e <9f2 x 

(2.7) /i(t,x) T / w (v)vn(x)rfv = -( /(v)v • n(x)dv. 

Jv-n(x)>0 Jvn(x)<0 

This boundary condition (2.4) guarantees the global conservation of mass [5]. 

In this paper we only apply a second order finite difference method to discretize the trans- 
port term of (2.1) but higher order schemes [13] may be applied. Then to keep the order 
of accuracy of the method, two ghost points should be added in each spatial direction. To 
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impose / at the ghost points, we will apply the inverse Lax-Wendroff procedure proposed 
in [21] for conservation laws. 

Suppose that the distribution function / at time level t n for all interior points are already 
known, we now construct / at the ghost points. 

2.1. One-dimensional case in space. We start with spatially one-dimensional problem, 
that is d x = 1. In this case the Boltzmann equation reads: 



x 



t 3 



(2-8) ^ + ^|^ = Iq( / ) ) ( x , v )e[x l: x r ] 

at ox e 

where xi and x r are the left and right boundaries respectively, v x is the component of phase 
space corresponding to x-direction. For the boundary condition in spatially one-dimensional 
case, the inward normal on the boundary in (2.4) is 



n(xi) = 



(1\ 



w 




To implement the numerical method, we assume the computational domain is a limited 
domain [a?min>^max] x [-V, V] 3 , where (xi,x r ) c [x m ; n , £ max ]. The computational domain is 
covered by a uniform Cartesian mesh x V/j, 



(2.9) 



— {x m i n — Xq < ••• < X{ < ••• < X nx — X max } , 

= {vj = j Av, j = (ji,...,j 3 ) eZ 3 , \j\<n v }. 



with the mesh size Ax and Av for space and velocity respectively. We only consider numerical 
method of ghost points near the left hand side boundary, since the procedure for right hand 
side boundary is the same. Figure 1 illustrates a portion of mesh near left boundary x\, which 
is located between xq and x\. 



-EEh 



X-l 



X 



x 2 



x 3 



S- 



v, 



3* 



Figure 1. A portion of mesh in spatially one dimensional case. • is interior 
point, ■ is ghost point, □ is the left hand side boundary. 



We construct / at each ghost point in following three steps: we perform an extrapolation of 
/ to compute a high order approximation of the outflow. Then, we compute an approximation 
of the distribution function at the boundary using Maxwell's boundary conditions. Finally, 
we apply the inverse Lax-Wendroff procedure for the inflow. 
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2.1.1. First step: Extrapolation of f for the outflow. At time t = t n we consider the outflow 
near the point xi, that is f(t,xi,\j) where Vj 1 < 0. We denote by fij an approximation of / 
at (x»,Vj) 

A natural idea is to extrapolate / at the left boundary X[ or the ghost points xq and x~\ 
using the values of / on interior points. For example from the values fij, f2j and faj, we 
can construct a Lagrange polynomial P2{x) e P2(M). Then by injecting xi, xq or X-\ into 
P2(x), we obtain the approximations of / at the ghost points and left boundary, i.e. fij, foj 
and f-ij- However, when a shock goes out of the boundary, the high order extrapolation 
may lead to a severe oscillation near the shock. To prevent this, we would like have a lower 
order accurate but more robust extrapolation. Therefore, a WENO type extrapolation [21] 
will be applied and described below (see subsection 2.3) for this purpose. 

2.1.2. Second step: Compute boundary conditions at the boundary. In the previous step, the 
outflow at the boundary is obtained by extrapolation. To compute the vlaues of / at the 
inflow boundary, we apply the Maxwell's boundary condition (2.4), i.e. 

(2.10) fij = (1-a) + aMtfij]. 

On the one hand the specular reflection portion is given straightly by the outflow at the left 
boundary, which is 

K[flj] = fij*, where 3* = (-31,32, ■ ■ ■ ,h)- 
On the other hand the diffuse one is computed by a half Maxwellian 

/ | v f 

M[fij] = mexpl—^r 
where T\ is the given temperature at the left wall and \x\ is given by 

m X v i - n ( x i) exp f _ fe^) = ~ £ v j - n ( x dfij- 

Vj-n(xi)>0 \ 1 I Vj-n(cc;)<0 

2.1.3. Third step: Approximation of f at the inflow boundary. Finally we compute the values 
of / at the ghost points for the inflow boundary. Here we cannot approximate / by an 
extrapolation, since the distribution function at interior points cannot predict the inflow. 
Then we extend the inverse Lax-Wendroff type procedure recently proposed in [11, 21, 23] 
for solving kinetic equations. At the left boundary x\, a first order Taylor expansion gives 

df 
dx 

Hence a first order approximation of / at ghost points is 

, s = -l,0. 

x=x t 

We already have fij in the second step, thus it remains to obtain an approximation of the 
first derivative. By reformulating (2.8), we have 

(2.12) §^ 

ox 



fj( x ) = fij + (x-xi) 



df 

(2.11) fsj - fij + (xs-xi) ~~j 



+ 0(Ax 2 ). 

X=Xl 



= -(-!-W)) 

v x \ at e 1 



X = Xl 



Now instead of approximating the first derivative d x f\ x=Xl , we compute the time derivative 
dtf\x=x t and the collision operator Q(/)|a;=a;i> An approximation of the time derivative can be 
computed by using several fij at previous time levels. Different approximation are obtained 
either a first order approximation reads 

fn _ rn—1 
Jlj Jlj 

At ' 



dt 



X=Xl 
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where At is the time step, or one can use a WENO type extrapolation to approximate the 
time derivative (see subsection 2.3 below). 
The last term Q(/)| 

x=xi can be computed explicitly by using fi j obtained in previous two 
steps. Clearly this procedure is independent of the values of / at interior points. 

Remark 2.1. Let us observe that when a = we have a pure specular reflection boundary 
condition. A mirror procedure can be used to approximate f at the ghost points. More 
precisely, by considering the boundary as a mirror, we approximate the distribution at the 
ghost points f(x s ,Vj) as 

f(x s ,Vj) = f(2xi-x s ,Vj*), where j* = {-jx,j2, ■ ■ ■ ,h) 

where 2 x\ - x s is the mirror image point of x s . Since 2 located in interior domain, 

we can approximate f(2xi -x s ,Vj+) by an interpolation procedure. 

2.2. Two-dimensional case in space. The previous approach can be generalized to spa- 
tially two-dimensional problem. We assume d x = 2 in equation (2.1) 

dl +Vx d£ +v d£ = l l 
dt x dx v dy e 

where the distribution function /(t,x, v) is defined in (i,x, v) e R + x x R 3 with x = (x,y). 
We consider a computational domain [a; m in, i m ax] x [ymin,2/max] x [^^j^] 3 , such that Q, c 
bmin^max] x [y min ,y max ] and /(t,x,v) 5? 0, for all ||v|| > V. 

The computational domain is covered by an uniform Cartesian mesh x V^, where X^, 
are defined similarly to (2.9). The mesh steps are respectively Ax, Ay and Av. In 
Figure 2, we present a portion of spatial mesh near the boundary. From a ghost point x g , we 
can find an inward normal n, which crosses the boundary at x p . 



(2-13) ^" + ^+^ = 7 G(/). 




Figure 2. Spatially two-dimensional Cartesian mesh. • is interior point, ■ is 
ghost point, □ is the point at the boundary, O is the point for extrapolation, 
the dashed line is the boundary. 



For the 2D case in space, the numerical approximation of the distribution function / at 
ghost points is similar to the one dimensional case. However, there are two major differences. 
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First to compute TZ[f] in the second step, the corresponding outflow may not locate on phase 
space mesh. Secondly to approximate the normal derivative in the third step, besides the 
time derivative and collision operator we need also the tangential derivative at x p . Once 
again, we present the method in three steps: 

2.2.1. First step: Extrapolation of f for outflow. Let us assume that the values of the distri- 
bution function / on the grid points in Q are given. To approximate / at a ghost point, for 
instance x s , we first construct a stencil £ composed of grid points of f2 for the extrapolation. 
For instance as it is shown in Figure 2, the inward normal n intersects the grid lines y - Ui y , 
yi y +i, yi y +2- Then we choose the three nearest point of the cross point in each line, i.e. marked 
by a large circle. From these nine points, we can build a Lagrange polynomial q2 (x) e Q2(M 2 ). 
Therefore we evaluate the polynomial ^(x) at x s or x p , and obtain an approximation of / 
at the boundary and at ghost points. As for the ID case, a WENO type extrapolation can 
be used to prevent spurious oscillations, which will be detailed in subsection 2.3. 

2.2.2. Second step: Compute boundary conditions at the boundary. In the previous step, we 
have obtained the outflow /(x p , v • n < 0) at the boundary x p . By using (2.4) as we did for 
the ID case, we can similarly compute the distribution function / for v • n > 0. However this 
time to compute the distribution function for specular reflection 

ft[/(x p , v)] = /(x p , v - 2(v ■ n)n), Vv e V fc) 

the vector fields v - 2 (v • n) n may not be located on a grid point. Therefore, we interpolate 
/ in phase space (x p , v-2(v-n)n) using the values computed from the outflow /(x p , v) such 
that v • n > 0. 

2.2.3. Third step: Approximation of f at the inflow boundary. We have obtained the values 
of / at the boundary points x p for all v e in previous two steps. Now we reconstruct the 
values of / for the velocity grid points such that v • n < at the ghost point x g by a simple 
Taylor expansion in the inward normal direction. To this end, we set up a local coordinate 
system at x p by 



x 



(x\ _ I cos 9 sm9\(x\ 
yj l-sin0 cos9J\yJ , 



where 9 is the angle between the inward normal n and the x-axis illustrated in Figure 2. 
Thus the first order approximation of /(x s , v) reads 

df 

/(x 9 ,v) a /(x p ,v) + (xg-Xp) — (£p,v), 

where /(x p ,v) = /(x p ,v) and ||(x p , v) is the first order normal derivative at the boundary 

x p . To approximate ||(x p ,v), we use inverse Lax-Wendroff procedure. Firstly, we rewrite 
the equation (2.13) in the local coordinate system as 

where v x = v x cos6> + v y sin9, v y = -v x sin9 + v y cos9. Then a reformulation of (2.14) yields 

I<-^-^(f-|->)L; 

Finally instead of approximating ||(x p ,v) directly, we approximate the time derivative |£, 
tangential derivative || and collision operator Q(f). Similarly as in spatially ID case, we 
compute |£ and Q(f). It remains to approximate ||. For this, some neighbor points of 
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x p at the boundary are required (See the empty squares in Figure 2). We then perform an 
essentially non-oscillatory (ENO) procedure [10] for this numerical differentiation to avoid 
the discontinuity. 



2.3. WENO type extrapolation. A WENO type extrapolation [21] was developed to pre- 
vent oscillations and maintain accuracy. The key point of WENO type extrapolation is to 
define smoothness indicators, which is designed to help us choose automatically between the 
high order accuracy and the low order but more robust extrapolation. Here we describe this 
method in spatially ID and 2D cases. Moreover we will give a slightly modified version of 
the method such that the smoothness indicators are invariant with respect to the scaling of 
/■ 



2.3.1. One- dimensional WENO type extrapolation. Assume that we have a stencil of three 
points £ = {xi,X2,xs} showed in Figure 1 and denote the corresponding distribution function 
by /i, /2, /3- Instead of extrapolating / at ghost point x g by Lagrange polynomial, we use 
following Taylor expansion 



_ 2 (x g -xi) 2 d k f 
Jg ~ 2-i 

k=0 



k\ 



dx k 



X=Xl 



We aim to obtain a (3 - k)-th order approximation of denoted by /j , k = 0, 1,2. 

Three candidate substencils are given by 

S r — {xi, . . . , x r+ i}, t — 0,1,2. 

In each substencil S r , we could construct a Lagrange polynomial p r (x) e P r (R) 

' PoO) = fi, 

PiO) = fi + — : — {x-xi), 
Ax 

P2{x) = fi + —r (x -xi) + — — {x - xi)(x - x 2 ). 

Ax 2Ax z 

We now look for the WENO type extrapolation in the form 



r=0 



where w r are the nonlinear weights depending on /j. We expect that has (3 - /c)-order 
accurate in the case f(x) is smooth in The nonlinear weights are given by 



Qy 



W r 



v^2 

L s =0&s 



with 



OLr - 



(e + /3,) 2 ' 
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where e = 10" 6 and f3 r are the new smoothness indicators determined by 



A) 



Ax 2 , 



+ fi+fiti Jx ° \ dx ) 



dx 



dx 



(/2-/1) 2 
£ + /? + /* 

1 

= -^7— 72T72-72t( 61 /i + 16 °/2 2 + 25/| + 74/x/ 3 - 192/x/ 2 - 124/ 2 / 3 ). 
We remark that the smoothness indicators (5\ and /3 2 have the factors — r l +1 2 , which 



ft = 



guarantee that the indicators are invariant of the scaling of fa 



2.3.2. Two- dimensional extrapolation. The two-dimensional extrapolation is a straightfor- 
ward expansion of ID case. The substencils S r , r = 0,1,2 for extrapolation are chosen 
around the inward normal n such that we can construct Lagrange polynomial of degree r. 
For instance in Figure 2, the three substencils are respectively 

' S = {(x ix ,y iy )}, 

51 = {{xi x -i,yi y ),(xi x ,yi y ),(xi x ,yi y+ i),(xi x+ i,yi y+ i)}, 

52 = {( x i x -l,yiy),( x i x ,yi y )i( x i x+ i,yi y ),(Xi x -i,yi y+ i), 
{ x ixiyi y + l)i { x ix + liyiy + l)i ( x i x iyiy+2)i { x ix+li yiy+2), ( x i x +2,yi y +2)}- 

Once the substencils S r are chosen, we could easily construct the Lagrange polynomials in 

Q r (R 2 ) 

r r 

9r(x) = Y, T, a l,rn xl y m 
m=0l=0 

satisfying 

q r (x) = /(x), x e S r . 
Then the WENO extrapolation has the form 

2 

(2.16) /(x) = Y, w r q r (x), x e S r , 

r=0 

where w r are the nonlinear weights, which are chosen to be 



with 



w r 



a r 



d r 

{e + f} r )< 



where e = 10 6 , do = Ax 2 + Ay 2 , d\ = \J Ax 2 + Ay 2 , d 2 = 1 - do _ di- Pr are the smoothness 
indicators determined by 

A) = Ax 2 + Ay 2 , 

Pr = : . ^ 1 — £ f \K\W- 1 (D a q r (x)) 2 dx,r = l,2, 



e + ExeSV /( x ) 2 l<\ a \<r ' 
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where a is a multi-index and K = [x p - Ax/2, x p + Ax/2] x [y p - Ay/2, y p + Ay/2], x p = (x p ,y p ). 



3. Application to the ES-BGK model 

The Boltzmann equation (2.1) governs well the evolution of density / in kinetic regime 
and also in the continuum regime [ ]. However the quadratic collision operator Q(f) has a 
rather complex form such that it is very difficult to compute. Hence different simpler models 
have been introduced. The simplest model is the so-called BGK model [ ], which is mainly 
a relaxation towards a Maxwellian equilibrium state 

(3.1) 2(/) = -CML7]-/), 

e 

where r depends on macroscopic quantities p and T. 

Although it describes the right hydrodynamical limit, the BGK model does not give the 
Navier-Stokes equation with correct transport coefficients in the Chapman-Enskog expansion. 
Holway et al. [ ] proposed the ES-BGK model, where the Maxwellian in the relaxation 

term of (3.1) is replaced by an anisotropic Gaussian G[f]- This model has correct conservation 
laws, yields the Navier-Stokes approximation via the Chapman-Enskog expansion with a 
Prandtl number less than one, and yet is endowed with the entropy condition [1]. In order to 
introduce the Gaussian model, we need further notations. Define the opposite of the stress 
tensor 

(3.2) e(t,x) = ± f (v-u)®(v-u)/(t,x,v)dv. 

p JR 3 

Therefore the translational temperature is related to the T = tr(0)/3. We finally introduce 
the corrected tensor 

T(t,x) = [(l-u)Tl + u@](t,x), 

which can be viewed as a linear combination of the initial stress tensor and of the isotropic 
stress tensor TI developed by a Maxwellian distribution, where I is the identity matrix. 

The ES-BGK model introduces a corrected BGK collision operator by replacing the local 
equilibrium Maxwellian by the Gaussian G[f] defined by 

G[f]= , P ej (v-u)T-*(v-u) 



y/det(2irT) 

Thus, the corresponding collision operator is now 

(3-3) Q(/) = T -{Q[f]-f), 

where r depends on p and T, the parameter -1/2 < v < 1 is used to modify the value of the 
Prandtl number through the formula 

2 1 

- < Pr = < +oo. 

3 1-1/ 

It follows from the above definitions that 



(3.4) 



X,v/(v)dv = / Rs vS[/](v) rf v = p„, 

I |2 I |2 

/ — f(v)dv= f ^-£7[/](v)cZv = E 



INVERSE LAX-WENDROFF METHOD FOR BOUNDARY CONDITIONS OF BOLTZMANN EQUATIONS 11 

and 

/ (v - u) <g> (v - u) /(v) dv = pQ, 

Jm 3 

* 

f (v-u)®(v-u)0[/]dv = pT. 
v Jm 3 

This implies that this collision operator does indeed conserve mass, momentum and energy 
as imposed. 

In this section, we will first recall the implicit-explicit (IMEX) scheme to the ES-BGK 
equation proposed in [3]. Then we apply our ILW procedure to treat the boundary condition 
for ES-BGK model case. 

3.1. An IMEX scheme to the ES-BGK equation. We now introduce the time dis- 
cretization for the ES-BGK equation (2.1), (3.3) 

df 



(3.5) 



+ v . Vx/ = _ (£[/] _ /)> xf!)c R d* v 6 M 3 
at e 



I /(0,x,v) = /o(x,v), xeO.vfR 3 , 

where r depends on p, u and T. 

The time discretization is an IMEX scheme. Since the convection term in (3.5) is not stiff, 
we will treat it explicitly. The source terms on the right hand side of (3.5) will be handled 
using an implicit solver. We simply apply a first order IMEX scheme, 



i f 



(3.6) 



n+1 



■n+1 



jn +v-V x /" = — (G[f n+1 ]-f n+1 ), 



This can be written as 

cn+1 



(3.7) 



/' 



At 

[ ,/°(x,v)=/ (x,v). 



r n+1 Ai 



■G[f n+l ], 



e + T n+1 At LJ ' J e + r n+1 At 

where G(f n+1 ) is the anisotropic Maxwellian distribution computed from f n+1 . Although 
(3.7) appears nonlinearly implicit, since the computation of / n+1 requires the knowledge of 
G[f n+1 ], it can be solved explicitly. Specifically, upon multiplying (3.7) by </>(v) defined by 



<Kv) == i 



-4) 



and use the conservation properties of Q and the definition of G[f] in (2.3), we define the 
macroscopic quantity U by U := (p,pu,E) computed from / and get 

F c T n+1 At r 



rn+1 



or simply 
(3.8) 



e + t 



Thus C/ n+1 can be obtained explicitly. This gives p n+1 ,u n+1 and T n+1 . Unfortunately, it is 
not enough to define G[f n+1 ] for which we need p n+1 G n+1 . Therefore, we define the tensor 
S by 

(3.9) S n+1 := J 3 v ® v / n+1 dv = p n+1 (9 n+1 + u n+1 ® u n+1 ) 

and multiply the scheme (3.7) by v ® v. Using the fact that 

/ v® v^[/](v)dv = p (T +u®u), 

JR 3 
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and (3.9), we get that 



(3.10) 



^n+l 



e + (1 - v) r r 



(l-u)r n+L At 
e + (I - u) r n+1 At 



— — - At f vgvvVxfrfv) 



P 



„n+l (yn+lj + u »+l 0u »tlj 



Now G[f n+1 ] can be obtained explicitly from U n+l and S n+1 and then / n+1 from (3.7). 
Finally the scheme reads 



(3.11) 



U 



n+l 



^n+1 



£ <Kv)(r-Atv-V x / n )dv, 

^ (V_At f v®vvV x / n dv) 

e+(l-v)r n+1 At\ Jr* I 



{l-u)T n+l At 
e + (l-u)r n+1 At ' 



■■" 1 (T n+1 I + u n+1 ®u" +1 ) 



/ 



71+1 



e + r n+1 Ai 



[/"-Atv-Vx/"] + 



r n+1 At 
e + r n+1 At 



Q[f n+1 l 



The scheme (3.11) is an AP scheme for (3.6). On the one hand, although (3.6) is nonlinearly 
implicit, is can be solved explicitly. On the other hand, the scheme (3.11) preserves the correct 
asymptotic [ ], which means when holding the mesh size and time step fixed and letting 
the Knudsen number go to zero, the scheme becomes a suitable scheme for the limiting 
hydrodynamic models. 



3.2. Inverse Lax-WendrofF procedure for boundary conditions. We have described 
the numerical method for boundary condition to general kinetic equations in spatially ID 
and 2D case. To implement this method, it remains to replace the collision operator Q(f) 
in (2.12) or (2.15) by the ES-BGK operator (3.3). 

Assume that the approximation to the distribution function at the boundary /(x p ,v) is 
known for all v e W Then, the macroscopic quantities p, u and T at the boundary point x p 
can be obtained using (2.3) and (3.4). Therefore, substituting these macroscopic quantities 
in (3.2), we compute the stress tensor O at the boundary point x p , such that the corrected 
tensor 7~(x p ). Thus G[f] is computed for all points (x p , v), where v e V^. 



4. Numerical examples 

In this section, we present a large variety of test cases in ld x and 2d x in space and three 
dimensional in velocity space showing the effectiveness of our method to get an accurate 
solution of Boltzmann type equations set in a complex geometry with different boundary 
conditions. We first give an example on a flow generated by gradients of temperature, which 
has already been treated by DSMC or other various methods [5]. 

Finally, we present some numerical results in 2d x . 

4.1. Smooth solutions. We consider the ES-BGK equation (2.1)-(3.3) 

l!t +Vx %= ^ Q(/) " ^ e (-°- 5 ' - 5 )' v6lR3 > 
- f(t = 0) = / , 
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with an initial datum /o which is a perturbation of the constant state in space and a 
Maxwellian distribution function in velocity, that is, 

/o(x ' V) = ^) 6XP (~~2~) ' X 6 ( "°- 5 ' °' 5) ' V 6 ^ 

with a density po - 1 + 0.1 cos(2-7rx). We consider purely diffusive boundary conditions with 
a wall temperature T w = 1. The solution is expected to be smooth for short time and then 
may develop a discontinuity at the boundary, which may propagate in the physical domain. 

We perform several numerical simulations on a time interval [0,t en d] with t en( i = 1, a 
computational domain in space Jo = [-7r/6,7r/6] such that (-1/2,1/2) c I and a domain in 
velocity V = [-8,8] 3 . Then, we choose a grid in space for Iq constituted of n x = n points 
and a grid for the velocity space with n v - n points for each direction with respectively 
n = 32, 64,..., n = 512. Let us emphasize that the boundary points x = -1/2 and x = 1/2 are 
not exactly located on a grid point. Since we don't know an exact solution of the problem, 
we compute relative errors. More precisely, an estimation of the relative error in L 1 norm at 
time T is given by 

e2h = \\fh(T)-f2h(T)\\ L i, 
where fh represents the approximation computed from a mesh of size h = (Ax,Av). The 
numerical scheme is said to be k-th order if &2h ^ C||^|| > f° r an < \h\ « 1. 

In Table 1 we compute the order of convergence in L 1 norm of our numerical methods. We 
can clearly see the expected second order convergence. Moreover, we verify experimentally 
that our scheme is also second-order accurate at the boundary since the discontinuity occuring 
at v x = is perfectly located. 



n x x n v 


L 1 error 


Order 


L 1 error at the boundary 


Order 


32 2 


8.8833 10" 4 


X 


3.909 10- 3 


X 


U z 


2.5221 10" 4 


1.94 


5.832 10~ 3 


X 


128 2 


6.551110^ 


1.88 


2.341 10" 4 


4.1 


256^ 


1.7829 1(T & 


1.91 


5.81110^ 


2.01 


512^ 


4.4571 l(T b 


2 


1.57310^ 


1.89 



Table 1. Smooth solutions: Experimental order of convergence in L norm. 
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Figure 3. Smooth solutions: Experimental order of convergence in L 1 norm 
(1) in the physical domain (2) at the boundary. 
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4.2. Flow generated by a gradient of temperature. We consider the ES-BGK equation 
(2.1)-(3.3), 

' % + Vx % = \ QU) ' S£ (- 1 / 2 . 1 /2),»«R 3 1 

1 / \v\ 2 \ 

fit = 0, x. v ) = — — exp — — , 
. ' 2nT \ 2T / 

with Tq(x) = 1 and we assume purely diffusive boundary conditions on x = -1/2 and x - 1/2, 
which can be written as 

f(t,x,v) = p(t,x) f w (v), if(x,v x ) e {-1/2} xM + and(x,^) e {1/2} x R" 

where \x is given by (2.7). This problem has already been studied in [2 ] using DSMC for 
the Boltzmann equation or using deterministic approximation using a BGK model for the 
Boltzmann equation in [16, 5]. 

Here we apply our numerical scheme with the ES-BGK operator (3.3) and choose a com- 
putational domain in space I Q = [-7r/6,7r/6] such that (-1/2,1/2) c I and [-8,8] 3 for the 
velocity space with a number grid points n = 32 in each direction and the time step At = 0.001. 

The main issue here is to capture the correct steady state for which the pressure is a 
perturbation of a constant state with a Knudsen layer at the boundary [16, 24]. 

In Figure 4, we represent the stationary solution (obtained approximately at time t en d = 25 
for e - 0.1 up to t en d - 75 for e = 0.025) of the temperature and the pressure profile. The 
results are in a qualitative good agreement with those already obtained in [24] with DSMC. 
More precisely, the boundary layer (Knudsen layer) appears in the density and temperature as 
well as the pressure, but it is small for all the quantities. The magnitude in the dimensionless 
density, temperature, and pressure is of order of e and the thickness of the layer is, say O(e). 
In the density and temperature profiles, we cannot observe it unless we magnify the profile in 
the vicinity of the boundary (see the zoom in Figure 4). Instead, since the pressure is almost 
constant in the bulk of the gas, we can observe perfectly the boundary layer by magnifying 
the entire profile. Let us emphasize that, as it is shown in Figure 4 the Knudsen layer is a 
kinetic effect, which disappears in the fluid limit (e ->■ 0). 

These results provide strong evidence that the present treatment of boundary conditions 
using WENO extrapolation and inverse Lax-Wendroff method can be used to determine the 
state of a gas under highly non-equilibrium conditions. Using deterministic methods, we can 
investigate the behavior of gases for situations in which molecular diffusion is important e.g., 
thermal diffusion. 

Also let us mention that a quantitative comparison between our results (3d v with ES-BGK 
operator) and [2 I] (3d v Boltzmann with hard sphere potential) or [16] (3d v BGK) gives a 
very good agreement on the values of the Knudsen layer and the values of the pressure inside 
the domain. 

4.3. High-speed flow through a trapezoidal channel. In this section we deal with spa- 
tially two-dimensional ES-BGK model in a trapezoidal domain. We attempt to get some 
steady state as 

v x^-+VyT- = -Q{f), 

ox ay e 

where v e Q x and v e M 3 . Here we will reproduce a numerical test performed in [19] but with 
our ILW method. The computational domain is a trapezoid 

S7 X = {x = (x, y), < x < a, < y < b + xtan(a)} 

as shown in Figure 5 for the parameters 

a = 2.0, b = 0.4, a = arctan(0.2). 
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Figure 4. Flow generated by a gradient of temperature: (1) temperature (2) 
pressure for various Knudsen numbers e = 0.025, 0.05 and 0.1. 



Boundary conditions are denned separately for each of the four straight pieces 

dn x = r, u r b u r r u r t 

denoting the left, bottom, right and top parts of the boundary respectively. The bottom part 
represents the axis of symmetry, so we use specular reflection (2.5) there, i.e. 

/(x,v) = (/(x,v-2(vn(x))n(x)), x e T b , v y > 0. 

On the right part we are modeling outflow (particles are permanently absorbed), i.e. 

(4.1) /(x,v) = 0, xer r , v x <0. 
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Figure 5. Trapezoidal domain f2 x . 



On the left part there is an incoming flux of particles, i.e. 
(4.2) /(x,v) = / in (x,v) =iWi„(v), xeT h v x >0, 

with an inflow Maxwellian 

j. ir i x An / |v~ ^n| 2 

On the top part of the boundary, we consider a diffuse reflection (2.5) of particles, with a 
Maxwellian distribution function 

Mr (v\ = PTfTll- 

2T t 



M rt (v) =exp(-Sr|. 



In the numerical experiments we assume 

Pin = 1, T in = 1, T t = 1.05 
and consider the inflow velocity in the form 

/1\ 

V in = Machjn^Tin , 

w 

where Mach; n = 5 and 7 = 1.4. 

To start the calculation we take an uniform initial solution equal to the values defined by 
the left boundary condition: 

fo(x,v) = - — /9 ' n ; - exp I - ■ V J^ 11 ^ I , x e $7 X , vfl 3 . 
J V ' ; (27rT in ) 3 / 2 F \ 2T in /' 

We define the Mach number from the macroscopic quantities, computing the moments of 
the distribution function with respect to v e R 3 , by 

Mach = 4^=, 

where c := V7T is the sound speed. 

We apply our inverse Lax-Wendroff method to the boundary conditions at 9f2 x . More 
precisely, we extrapolate first the outflow at ghost points corresponding to the four straight 
pieces. Then we impose directly the inflow at the boundaries Ti and T r by (4.1), (4.2), since 
they are independent of outflow. While the inflow of Tf, and Tj is computed by specular and 
diffuse reflection. Finally we use inverse Lax-Wendroff procedure to compute inflow at ghost 
points. 
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Figure 6. Density with e = 5 
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Figure 7. Temperature, e = 5 
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Figure 8. Mach with e = 5 



In following the sequel, numerical experiments are performed on a mesh of size 96 x 48 on 
space domain f] x . For velocity space we choose limit domain [-12,12] x [-8,8] x [-8,8] with 
the grid point number as 64 x 48 x 12. Moreover for the ES-BGK operator (3.3) we choose 
v = -0.5. We first consider the weak collision case, i.e. Kn = 5. In Figures 6-8, we show on 
the left hand side, the contour plots of the density, the temperature and the Mach number 
while the right hand side plots show the absolute values of these quantities plotted along 
the axis of symmetry y = 0. We observe that the flow changes when we consider different 
Knudsen numbers Kn = 0.05 and Kn = 5. The corresponding results are shown in Figures 9- 
11. The significant difference between these two case can be observed in Mach number. In 
the case Kn = 5, the Mach number reaches its maximum at T r while in the case Kn = 0.05 its 
maximum is at T[. We can observe also in the case Kn = 0.05 that there is a clear maximum 
of the density in the middle of the domain. In the same region the temperature reaches its 
maximum. 

4.4. High-speed flow around an object. In this section, we desire to simulate high-speed 
airflow around a half airfoil (see Figure 12). The boundary is separated by four parts 



do x = rjur fe ur r ur 4 . 
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Figure 9. Density with e = 0.05 



1.25 
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Figure 10. Temperature with e = 0.05 




Figure 11. Mach with e = 0.05 



On the right and left hand sides T r , Ti, we use the same incoming flux (4.1), (4.2). On the 
top part Tt, the incoming flow is given by the initial value 

/(x,v)= Pin NQ/0 exp(- |V "^ n|2 V xeT(, v e M 3 . 
JK ' (2vrT in )3/2 p \ 2T in /' 

Finally at the bottom Fb, we use a purely specular reflection boundary condition. The 
parameters p- m , T; n , T t , 7, v have the same values as in the previous test. We use again (4.3) 
as the initial solution. 

We note that on the profile of airfoil we cannot use the neighbor points to approximate the 
tangential derivative ^£ in (2.15). It is because these neighbor points are not on the same 
straight. Here we approximate the tangential derivative by using the distribution function of 
interior domain. 

In the following tests, we consider only the situations in hydrodynamic regime, i.e. e = 0.05, 
for comparing the ones in literature [8, 14]. A 150 x 100 mesh is used in domain f2 x . We use 
a limit velocity domain [-8, 8] 3 with mesh size 48 x 48 x 12. Two different tests of transonic 
airflow around this half airfoil are considered: Machj n < 1 and Machi n > 1. 

We choose first Machi n = 0.85. So we observe in Figures 13-15 that the flow field around 
the object includes both sub- (Mach < 1) and supersonic (Mach > 1.2) parts. The transonic 
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Figure 12. Domain including a half airfoil ri x . 
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FIGURE 13. Density with Mach in = 0.85 

(0.8 < Mach < 1.2) period begins when first zones of Mach > 1 flow appear around the 
object. Supersonic flow can decelerate back to subsonic before the trailing edge. In the case 
Machj n = 1.2 (see in Figures 16-18), the zone of Mach > 1 flow increases towards both leading 
and trailing edges. There is a normal shock created at trailing edge. The flow decelerates 
over the shock, but remains supersonic. Moreover a normal shock is created ahead of the 
object, and the only subsonic zone in the flow field is a small area around the object's leading 
edge. 

5. Conclusion 

In this paper we present an accurate method based on Cartesian mesh to deal with complex 
geometry boundary for kinetic models set in a complex geometry. We desire to reconstruct 
the distribution function / on some ghost points for computing transport operator. For this 
we proceed in three steps: first we extrapolate the distribution function / on ghost points 
for outflow. Then we use the boundary conditions to compute the inflow at the boundary. 
Finally we implement an inverse Lax-Wendroff procedure to give an accurate approximation 



1 - 
0.5- 
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Figure 14. Temperature with Mach in = 0.85 
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Figure 15. Mach with Mach in = 0.85 

of / for inflow on the ghost points. A spatially one-dimensional example is given to show that 
this method has second order accuracy in L 1 norm. Moreover some ID x 3D and 2D x 3D 
illustrate that our method can reproduce the similar results as the ones in literature. 
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